/*
 heat-spectral lib. solves time-dependent PDEs using SEM-NI method.
 Copyright (C) 2009, Altieri Michele and Stefanoni Andrea
 
 This program is free software: you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation, either version 3 of the License, or
 (at your option) any later version.
 
 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 GNU General Public License for more details.
 
 You should have received a copy of the GNU General Public License
 along with this program. If not, see <http://www.gnu.org/licenses/>.
*/

#ifndef TUBE_HPP
#define TUBE_HPP

#include<cmath>
#include<spectral>

// *******************************************************************
// Diffusion-Advection Problem
// *******************************************************************
// u_t - \mu \Delta u + b\cdot\Nabla u = f(t)	in \Omega \times [0,T]
//
// u(x,0) = u0					in \Omega
// u(x,y) = uD					on \Gamma_D
// \Nabla u(x,t) \cdot n = 0			on \Gamma_N
// *******************************************************************

#define TUBE_PROBLEM TubeBilinear<spectral::Mesh>, Functional<spectral::Mesh,TubeFunc,true>

// Parameters of the physical problem
double length(8);
double width(2);
double tube(1);
double thickness1(0.1);
double thickness2(0.2);
double thickness = thickness1 + thickness2;
double radius = tube-thickness;

// Acceleration parameter
double speed(3600);

struct TubeFunc {
    double operator()(double x,double y,double t) const {
        if( (std::abs(y) > radius)
		&& (std::abs(y) < (tube-thickness2))
		&& (((x > 0.2*length) && (x < 0.4*length)) ||
		    ((x > 0.6*length) && (x < 0.8*length)))
		&& (t < 48) ) {
            return 20*std::sin(2.*M_PI/24*t);
        } else {
            return 0;
        }
    }
};

struct Tubemu {
    double operator()(double x, double y) const {
        if (std::abs(y) > tube) {
            return 2.216e-5*speed;	// Air
        } else if (std::abs(y) > (radius+thickness1)){
            return 2.8e-7*speed;	// Glass wool
        } else if (std::abs(y) > radius){
            return 1.12e-4*speed;	// Copper
        } else {
            return 2.204e-5*speed;	// Nitrogen
        }
    }
};

struct Tubeb1 {
    double operator()(double x, double y) const {
	if (std::abs(y) <= radius) {
	    return 0.1*(radius-y)*(radius+y)*speed;
	} else {
	    return 0;
	}
    }
};

struct Tubeb2 {
    double operator()(double x, double y) const {
	return 0;
    }
};

template<class MESH>
class TubeBilinear {
    public:
	TubeBilinear(MESH const & me): me(me) {
	    for(int i=0;i<4;++i) {
		bcond[i]=spectral::dirichlet;
	    }
	    bcond[1]=spectral::neumann;
	}

	template<class U,class V> inline
	    double operator()(U const & u,V const & v,int n) const {
		return
		    + me.prodL2(me.restriction(mu,n),me.Dx(u,n),me.Dx(v,n),n)
		    + me.prodL2(me.restriction(mu,n),me.Dy(u,n),me.Dy(v,n),n)
		    + me.prodL2(me.restriction(b1,n),me.Dx(u,n),v,n)
		    + me.prodL2(me.restriction(b2,n),me.Dy(u,n),v,n);
	    }

	// Initial data
	struct initial{
	    double operator()(double const & x, double const & y) const {
		return 300;
	    }
	};

	// Boundary conditions
	spectral::bctype bcond[4];

	// Function for the Dirichlet boundary condition
	struct dirichlet{
	    double operator()(double const & x, double const & y) const {
		return 300;
	    }
	};

	// Function for the Neumann boundary condition
	struct neumann{
	    double operator()(double const & x, double const & y) const {
		return 0;
	    }
	};

	// Methods on
	template<unsigned int N, typename T=void>
	struct on : public dirichlet {};
	
	template<typename T>
	struct on<1,T> : public neumann {};

    private:
	MESH const & me;
	Tubemu mu;
	Tubeb1 b1;
	Tubeb2 b2;
};

const spectral::Point<double> LA(0.,-width);
const spectral::Point<double> LB(length,width);

#endif
